Interpreting and Visualizing GAMs

knitr::opts_chunk$set(echo = TRUE,dpi = 400)

Significance and linearity

library(gamair)
data("mpg", package="gamair")

library(mgcv)

## Loading required package: nlme

## Warning: package 'nlme' was built under R version 3.6.2

## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

# Fit the model
mod_city4 <- gam(city.mpg ~ s(weight) + s(length) + s(price) + s(rpm) + s(width),
                 data = mpg, method = "REML")

# View the summary
summary(mod_city4)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## city.mpg ~ s(weight) + s(length) + s(price) + s(rpm) + s(width)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   25.201      0.188     134   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df      F  p-value    
## s(weight) 5.620  6.799 17.524  < 2e-16 ***
## s(length) 2.943  3.759  0.904    0.420    
## s(price)  1.000  1.000 16.647 6.68e-05 ***
## s(rpm)    7.751  8.499 16.486  < 2e-16 ***
## s(width)  1.003  1.005  0.006    0.939    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.831   Deviance explained = 84.7%
## -REML = 496.47  Scale est. = 7.0365    n = 199

Plotting the motorcycle crash data and model

mcycle <- MASS::mcycle

# Fit the model
mod <- gam(accel ~ s(times), data = mcycle, method = "REML")

# Make the plot with residuals
plot(mod,residuals=T)

# Change shape of residuals
plot(mod, residuals=T,pch=1,cex=1)

Plotting multiple auto performance variables

data("mpg", package="gamair")

# Fit the model
mod <- gam(hw.mpg ~ s(weight) + s(rpm) + s(price) + comp.ratio, 
           data = mpg, method = "REML")

# Plot the price effect
plot(mod,select=c(3))

# Plot all effects
plot(mod,pages=1,all.terms=T)

Visualizing auto performance uncertainty

# Fit the model
mod <- gam(hw.mpg ~ s(weight) + s(rpm) + s(price) + comp.ratio, 
           data = mpg, method = "REML")

# Plot the weight effect with colored shading
plot(mod, select = 1, shade=T, shade.col="hotpink")

# Make another plot adding the intercept value and uncertainty
plot(mod, select = 1, shift=coef(mod)[1])

Reading model diagnostics

set.seed(0)
dat <- gamSim(1,n=200)

## Gu & Wahba 4 term additive model

library(mgcv)
# Fit the model
mod <- gam(y ~ s(x0, k = 5) + s(x1, k = 5) + s(x2, k = 5) + s(x3, k = 5),
           data = dat, method = "REML")

# Run the check function
gam.check(mod)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-0.0001190691,0.0001259251]
## (score 460.9549 & scale 5.229925).
## Hessian positive definite, eigenvalue range [0.0001190726,97.53256].
## Model rank =  17 / 17 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##         k'  edf k-index p-value   
## s(x0) 4.00 2.54    0.99    0.47   
## s(x1) 4.00 2.25    0.95    0.18   
## s(x2) 4.00 3.94    0.84    0.01 **
## s(x3) 4.00 1.00    0.99    0.40   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fixing problems with model diagnostics

set.seed(0)
dat <- mgcv::gamSim(1,n=600, scale=0.6, verbose=FALSE)

# Fit the model
mod <- gam(y ~ s(x0, k = 3) + s(x1, k = 3) + s(x2, k = 3) + s(x3, k = 3),
           data = dat, method = "REML")

# Check the diagnostics
gam.check(mod)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-2.133524e-06,2.107185e-06]
## (score 1262.659 & scale 3.777109).
## Hessian positive definite, eigenvalue range [0.2592108,297.5029].
## Model rank =  9 / 9 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##         k'  edf k-index p-value    
## s(x0) 2.00 1.99    0.98    0.32    
## s(x1) 2.00 1.95    1.07    0.94    
## s(x2) 2.00 2.00    0.11  <2e-16 ***
## s(x3) 2.00 1.72    1.03    0.78    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Refit to fix issues
mod2 <- gam(y ~ s(x0, k = 3) + s(x1, k = 3) + s(x2, k = 12) + s(x3, k = 3),
           data = dat, method = "REML")

# Check the new model
gam.check(mod2)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-3.827e-08,3.651907e-08]
## (score 615.861 & scale 0.3924265).
## Hessian positive definite, eigenvalue range [0.3995874,297.5842].
## Model rank =  18 / 18 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##          k'   edf k-index p-value
## s(x0)  2.00  2.00    1.07    0.96
## s(x1)  2.00  2.00    0.97    0.20
## s(x2) 11.00 10.79    1.01    0.56
## s(x3)  2.00  1.89    0.98    0.29

Examining the overall concurivity in auuuto data

set.seed(0)

# Fit the model
mod <- gam(hw.mpg ~ s(length) + s(width) + s(height) + s(weight),
           data = mpg, method = "REML")

# Check overall concurvity
concurvity(mod,full=T)

##                  para s(length)  s(width) s(height) s(weight)
## worst    1.079374e-20 0.9303404 0.9322887 0.6723705 0.9603887
## observed 1.079374e-20 0.7534619 0.8757513 0.4869308 0.8793300
## estimate 1.079374e-20 0.8353324 0.7943374 0.4452676 0.8567519

Examining concurvity between auto variables

# Fit the model
mod <- gam(hw.mpg ~ s(length) + s(width) + s(height) + s(weight),
           data = mpg, method = "REML")

# Check pairwise concurvity
concurvity(mod,full=F)

## $worst
##                   para    s(length)     s(width)    s(height)    s(weight)
## para      1.000000e+00 4.799804e-26 5.458174e-21 4.926340e-23 3.221614e-25
## s(length) 4.759962e-26 1.000000e+00 8.336513e-01 6.058015e-01 8.797217e-01
## s(width)  5.458344e-21 8.336513e-01 1.000000e+00 4.099837e-01 8.953662e-01
## s(height) 4.927251e-23 6.058015e-01 4.099837e-01 1.000000e+00 3.665831e-01
## s(weight) 3.233688e-25 8.797217e-01 8.953662e-01 3.665831e-01 1.000000e+00
## 
## $observed
##                   para    s(length)     s(width)    s(height)    s(weight)
## para      1.000000e+00 1.128295e-29 4.467995e-32 9.887661e-34 6.730965e-31
## s(length) 4.759962e-26 1.000000e+00 7.511142e-01 2.827977e-01 8.232449e-01
## s(width)  5.458344e-21 5.077384e-01 1.000000e+00 1.186126e-01 7.813743e-01
## s(height) 4.927251e-23 2.284116e-01 3.313152e-01 1.000000e+00 2.900361e-01
## s(weight) 3.233688e-25 6.052819e-01 7.863555e-01 1.494913e-01 1.000000e+00
## 
## $estimate
##                   para    s(length)     s(width)    s(height)    s(weight)
## para      1.000000e+00 1.564968e-28 1.740649e-23 3.448567e-25 1.481483e-27
## s(length) 4.759962e-26 1.000000e+00 6.415191e-01 2.271285e-01 7.209033e-01
## s(width)  5.458344e-21 6.477497e-01 1.000000e+00 1.054762e-01 7.241891e-01
## s(height) 4.927251e-23 3.303484e-01 2.644827e-01 1.000000e+00 2.669300e-01
## s(weight) 3.233688e-25 7.235198e-01 6.913221e-01 1.390568e-01 1.000000e+00